Particle-COMETS: From dynamic FBA to spatial simulations with particles¶
Eran Agmon, University of Connecticut
This notebook is dedicated to developing a particle-COMETS simulation step by step using Vivarium. Each section of this notebook incrementally adds layers to our model, illustrating the different computational methods and how they integrate.
In [1]:
import numpy as np
# spatio-flux's customized Vivarium
from spatio_flux import SpatioFluxVivarium, render_path
Make a Vivarium¶
view available types and processes
In [2]:
# make a fresh vivarium
vi = SpatioFluxVivarium()
In [3]:
# view the available types
vi.get_types()
Out[3]:
['tuple', 'array', '', 'length^1_5*mass^0_5/time', 'length/time', 'enum', 'temperature', 'mass^0_5/length^1_5', 'length^2*mass/time^3', 'mass/length^3', 'length*mass/time^2', 'mass/length*time^2', 'length^2*mass/current^2*time^2', 'current*time', 'edge', 'map', 'interval', 'printing_unit/length', 'particle', 'mass/time^2', 'length*time/mass', 'wires', 'luminosity', 'length^2/time^2', 'emitter_mode', 'current^2*time^4/length^2*mass', 'current^2*time^3/length^2*mass', 'length^3/mass*time^2', 'mass/length*time', 'length^3*mass/current^2*time^4', 'substance/length^3', 'composite', 'union', 'quote', 'length^2*mass/current^2*time^3', 'length', 'length^1_5*mass^0_5/time^2', 'length^4*mass/time^3', 'current^2*time^4/length^3*mass', 'length^2*mass/current*time^3', 'time^2/length', 'schema', 'substrate_role', 'current*time/substance', 'length^2', 'number', 'length^2*mass/current*time^2', 'float', 'printing_unit', 'length/time^2', 'length^3/time', 'boolean', 'string', 'length^2/time', 'length^4*mass/current*time^3', 'boundary_side', 'length^2*mass/time', 'mass^0_5/length^0_5*time', '/time', 'bounds', '/temperature*time', 'maybe', 'tree', 'length*temperature', 'length^0_5*mass^0_5/time', 'length^0_5*mass^0_5', 'substance/time', 'process', 'mass/length', 'protocol', 'length^2*mass/time^2', 'positive_float', 'current*length^2', '/substance', 'mass/temperature^4*time^3', 'time/length', 'time', 'current*length*time', 'current*time/mass', 'list', 'mass', 'step', 'reaction', 'length^2*mass/substance*temperature*time^2', 'length/mass', 'integer', 'length*mass/current*time^3', 'path', '/printing_unit', '/length', 'substance', 'length^2*mass/temperature*time^2', 'current', 'kinetics', 'current*time^2/length^2*mass', 'any', 'mass/time^3', 'current*length^2*time', 'luminosity/length^2', 'length^3', 'mass/current*time^2', 'length*mass/current^2*time^2']
In [4]:
# view the available processes
vi.get_processes()
Out[4]:
['composite', 'ram-emitter', 'DynamicFBA', 'json-emitter', 'MinimalParticle', 'console-emitter', 'Particles', 'DiffusionAdvection']
dFBA¶
Dynamic Flux Balance Analysis (dFBA) extends traditional Flux Balance Analysis (FBA) to model the dynamic behavior of metabolic networks over time, allowing for the simulation of growth and substrate consumption in a changing environment.
In [5]:
# inspect the DynamicFBA process config
vi.process_config('DynamicFBA')
Out[5]:
{'model_file': 'string',
'kinetic_params': 'map[tuple[float,float]]',
'substrate_update_reactions': 'map[string]',
'biomass_identifier': 'string',
'bounds': 'map[bounds]'}
In [8]:
# dfba_config = vi.process_config('DynamicFBA', dataclass=True) # TODO get dataclass to configure
dfba_config = {
"model_file": "textbook",
"kinetic_params": {
"glucose": (0.5, 1),
"acetate": (0.5, 2)},
"substrate_update_reactions": {
"glucose": "EX_glc__D_e",
"acetate": "EX_ac_e"},
"biomass_identifier": "biomass",
"bounds": {
"EX_o2_e": {"lower": -2, "upper": None},
"ATPM": {"lower": 1, "upper": 1}}}
In [9]:
# make a fresh vivarium
v1 = SpatioFluxVivarium()
# add a dFBA process
v1.add_process(
name="dFBA",
process_id="DynamicFBA",
config=dfba_config)
v1.diagram(dpi='70')
Out[9]:
In [10]:
mol_ids = ["glucose", "acetate", "biomass"]
# add the molecular fields
for mol_id in mol_ids:
v1.add_object(
name=mol_id,
path=['fields'],
value=np.random.rand())
v1.connect_process(
name="dFBA",
inputs={"substrates": {
mol_id: ['fields', mol_id]
for mol_id in mol_ids}},
outputs={"substrates": {
mol_id: ['fields', mol_id]
for mol_id in mol_ids}})
# add an emitter to save results
v1.add_emitter()
v1.diagram(dpi='70', show_values=True)
Out[10]:
In [11]:
v1.set_value(path=['fields', 'glucose'], value=10)
v1.set_value(path=['fields', 'biomass'], value=0.1)
field = v1.get_value(['fields'])
print(field)
{'glucose': 10, 'acetate': 0.7464969100996626, 'biomass': 0.1}
In [12]:
# save a file with the exact simulation state
v1.save(filename='dFBA_t0')
Saved file: out/dFBA_t0.json
In [13]:
vx = SpatioFluxVivarium(document='out/dFBA_t0.json')
vx.run(interval=10)
vx.diagram(dpi='70', show_values=True)
Out[13]:
In [14]:
# run the simulation
v1.run(interval=60)
In [15]:
# view the timeseries
v1.get_timeseries(as_dataframe=True)
Out[15]:
| /global_time | /fields/glucose | /fields/acetate | /fields/biomass | |
|---|---|---|---|---|
| 0 | 0.0 | 10.000000 | 0.746497 | 0.100000 |
| 1 | 1.0 | 9.904762 | 0.758605 | 0.107943 |
| 2 | 2.0 | 9.802006 | 0.771581 | 0.116515 |
| 3 | 3.0 | 9.691146 | 0.785476 | 0.125766 |
| 4 | 4.0 | 9.571550 | 0.800341 | 0.135747 |
| ... | ... | ... | ... | ... |
| 56 | 56.0 | 0.000000 | 0.000000 | 0.986275 |
| 57 | 57.0 | 0.000000 | 0.000000 | 0.986275 |
| 58 | 58.0 | 0.000000 | 0.000000 | 0.986275 |
| 59 | 59.0 | 0.000000 | 0.000000 | 0.986275 |
| 60 | 60.0 | 0.000000 | 0.000000 | 0.986275 |
61 rows × 4 columns
In [16]:
# plot the timeseries
v1.plot_timeseries(
subplot_size=(8, 3),
combined_vars=[
[ # combine the variables into a single subplot
'/fields/glucose',
'/fields/acetate',
'/fields/biomass'
]])
Out[16]:
Spatial dFBA¶
In [17]:
mol_ids = ["glucose", "acetate", "biomass"]
rows = 2
columns = 1
# make a fresh vivarium
v2 = SpatioFluxVivarium()
for mol_id in mol_ids:
v2.add_object(
name=mol_id,
path=['fields'],
value=np.random.rand(rows, columns))
# add a dynamic FBA process at every location
for i in range(rows):
for j in range(columns):
dfba_name = f"dFBA[{i},{j}]"
v2.add_process(
name=dfba_name,
process_id="DynamicFBA",
config=dfba_config)
v2.connect_process(
name=dfba_name,
inputs={"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids}},
outputs={"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids}})
# add an emitter to save results
v2.add_emitter()
v2.diagram(dpi='70')
Out[17]:
In [18]:
# change some initial values
v2.merge_value(path=['fields', 'glucose', 0, 0], value=10.0)
v2.merge_value(path=['fields', 'biomass', 0, 0], value=0.1)
field = v2.get_value(['fields'])
print(field)
{'glucose': array([[10. ],
[ 0.61243865]]), 'acetate': array([[0.32678828],
[0.96336482]]), 'biomass': array([[0.1 ],
[0.41108596]])}
In [19]:
# run a simulation
v2.run(60)
In [20]:
# get a list of all the paths so they can be plotted together in a single graph
all_paths = [
[render_path(['fields', mol_id, i, j]) for mol_id in mol_ids]
for i in range(rows)
for j in range(columns)]
# plot the timeseries
v2.plot_timeseries(
subplot_size=(8, 3),
combined_vars=all_paths)
Out[20]:
In [21]:
v2.plot_snapshots()
Diffusion/Advection¶
This approach models the physical processes of diffusion and advection in two dimensions, providing a way to simulate how substances spread and are transported across a spatial domain, essential for understanding patterns of concentration over time and space.
In [22]:
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
'glucose',
'acetate',
'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, -0.1)}
# make a fresh Vivarium
v3 = SpatioFluxVivarium()
# add fields for all the molecules
for mol_id in mol_ids:
v3.add_object(
name=mol_id,
path=['fields'],
value=np.random.rand(n_bins[0], n_bins[1]))
# add a spatial diffusion advection process
v3.add_process(
name='diffusion_advection',
process_id='DiffusionAdvection',
config={
'n_bins': n_bins,
'bounds': bounds,
'default_diffusion_rate': diffusion_rate,
'default_diffusion_dt': diffusion_dt,
'advection_coeffs': advection_coeffs},
inputs={'fields': ['fields']},
outputs={'fields': ['fields']})
# add an emitter to save results
v3.add_emitter()
v3.diagram(dpi='70')
Out[22]:
In [23]:
v3.run(60)
In [24]:
v3.show_video()
COMETS (Computation Of Microbial Ecosystems in Time and Space)¶
COMETS combines dynamic FBA with spatially resolved physical processes (like diffusion and advection) to simulate the growth, metabolism, and interaction of microbial communities within a structured two-dimensional environment, capturing both biological and physical complexities.
In [25]:
bounds = (20.0, 10.0) # Bounds of the environment
n_bins = (12, 6)
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, 0.1)}
# make a fresh vivarium
v4 = SpatioFluxVivarium()
# initialize the molecular fields
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc
acetate_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1 # place some biomass
v4.add_object(
name='glucose',
path=['fields'],
value=glc_field)
v4.add_object(
name='biomass',
path=['fields'],
value=biomass_field)
v4.add_object(
name='acetate',
path=['fields'],
value=acetate_field)
# add a diffusion/advection process
v4.add_process(
name='diffusion_advection',
process_id='DiffusionAdvection',
config={
'n_bins': n_bins,
'bounds': bounds,
'default_diffusion_rate': diffusion_rate,
'default_diffusion_dt': diffusion_dt,
'advection_coeffs': advection_coeffs },
inputs={'fields': ['fields']},
outputs={'fields': ['fields']})
# add a dynamic FBA process at every location
for i in range(n_bins[0]):
for j in range(n_bins[1]):
dfba_name = f"dFBA[{i},{j}]"
v4.add_process(
name=dfba_name,
process_id="DynamicFBA",
config=dfba_config )
v4.connect_process(
name=dfba_name,
inputs={"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids} },
outputs={"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids}})
# add an emitter to save results
v4.add_emitter()
v4.diagram(dpi='70',
remove_nodes=[f"/dFBA[{i},{j}]" for i in range(n_bins[0]-1) for j in range(n_bins[1])])
Out[25]:
In [26]:
v4.run(60)
In [27]:
v4.show_video()
In [28]:
v4.plot_timeseries(
subplot_size=(8, 3),
query=[
'/fields/glucose/0/0',
'/fields/acetate/0/0',
'/fields/biomass/0/0'],
combined_vars=[[
'/fields/glucose/0/0',
'/fields/acetate/0/0',
'/fields/biomass/0/0']
])
Out[28]:
Particles¶
In [29]:
import numpy as np
from spatio_flux import SpatioFluxVivarium
In [30]:
vi.process_config('Particles')
Out[30]:
{'bounds': 'tuple[float,float]',
'n_bins': 'tuple[integer,integer]',
'diffusion_rate': {'_type': 'float', '_default': 0.1},
'advection_rate': {'_type': 'tuple[float,float]', '_default': (0, 0)},
'add_probability': 'float',
'boundary_to_add': {'_type': 'list[boundary_side]',
'_default': ['left', 'right']},
'boundary_to_remove': {'_type': 'list[boundary_side]',
'_default': ['left', 'right', 'top', 'bottom']}}
In [31]:
bounds = (10.0, 20.0) # Bounds of the environment
n_bins = (20, 40) # Number of bins in the x and y directions
v5 = SpatioFluxVivarium()
v5.add_process(
name='particle_movement',
process_id='Particles',
config={
'n_bins': n_bins,
'bounds': bounds,
'diffusion_rate': 0.1,
'advection_rate': (0, -0.1),
'add_probability': 0.3,
'boundary_to_add': ['top']})
v5.connect_process(
name='particle_movement',
inputs={'fields': ['fields'],
'particles': ['particles']},
outputs={'fields': ['fields'],
'particles': ['particles']})
v5.initialize_process(
path='particle_movement',
config={'n_particles': 2})
v5.add_emitter()
v5.diagram(dpi='70')
Out[31]:
In [32]:
v5.save('particle_movement')
Saved file: out/particle_movement.json
In [33]:
v5.run(100)
v5_results = v5.get_results()
In [34]:
v5.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif
In [35]:
v5.diagram(dpi='70')
Out[35]:
In [36]:
v5.save(filename='v5_post_run.json', outdir='out')
Saved file: out/v5_post_run.json
Minimal particle process¶
In [37]:
import numpy as np
from spatio_flux import SpatioFluxVivarium,render_path
from spatio_flux.processes.particles import get_minimal_particle_composition
In [38]:
bounds = (10.0, 20.0) # Bounds of the environment
n_bins = (20, 40) # Number of bins in the x and y directions
v6 = SpatioFluxVivarium()
# make two fields
v6.add_object(name='glucose',
path=['fields'],
value=np.ones((n_bins[0], n_bins[1])))
v6.add_object(
name='acetate',
path=['fields'],
value=np.zeros((n_bins[0], n_bins[1])))
v6.add_process(
name='particle_movement',
process_id='Particles',
config={
'n_bins': n_bins,
'bounds': bounds,
'diffusion_rate': 0.1,
'advection_rate': (0, -0.1),
'add_probability': 0.3,
'boundary_to_add': ['top']})
v6.connect_process(
name='particle_movement',
inputs={
'fields': ['fields'],
'particles': ['particles']},
outputs={
'fields': ['fields'],
'particles': ['particles']})
# add a process into each particles schema
minimal_particle_config = {
'reactions': {
'grow': {
'glucose': {
'vmax': 0.01,
'kcat': 0.01,
'role': 'reactant'},
'acetate': {
'vmax': 0.001,
'kcat': 0.001,
'role': 'product'}
}}}
particle_schema = get_minimal_particle_composition(v6.core, minimal_particle_config)
v6.merge_schema(path=['particles'], schema=particle_schema['particles'])
# add particles to the initial state
v6.initialize_process(
path='particle_movement',
config={'n_particles': 1})
v6.diagram(dpi='70')
Out[38]:
In [40]:
# v6.set_value('particles', {})
# v6.initialize_process(
# path='particle_movement',
# config={'n_particles': 1})
v6.run(60)
v6_results = v6.get_results()
In [41]:
v6.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif
add diffusion¶
In [52]:
document = v6.make_document()
v6x = SpatioFluxVivarium(document=document)
particle_id = list(v6x.composite.state['particles'].keys())[0]
# v6x.composite.composition['particles'] = {}
# v6x.composite.composition['particles']
# v6x.composite.state['particles'][particle_id]['minimal_particle']['_type'] = 'process'
v6x.diagram(dpi='70')
Out[52]:
In [50]:
document = v6.make_document()
v6x = SpatioFluxVivarium(document=document)
In [51]:
v6x.run(60)
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) Cell In[51], line 1 ----> 1 v6x.run(60) File ~/code/vivarium-interface/vivarium/vivarium.py:574, in Vivarium.run(self, interval) 571 if not self.emitter_paths: 572 self.add_emitter() --> 574 self.composite.run(interval) File ~/code/process-bigraph/process_bigraph/composite.py:800, in Composite.run(self, interval, force_complete) 798 for path in self.process_paths: 799 process = get_path(self.state, path) --> 800 full_step = self.run_process( 801 path, 802 process, 803 end_time, 804 full_step, 805 force_complete) 807 # apply updates based on process times in self.front 808 if full_step == math.inf: 809 # no processes ran, jump to next process File ~/code/process-bigraph/process_bigraph/composite.py:697, in Composite.run_process(self, path, process, end_time, full_step, force_complete) 694 if state is None: 695 breakpoint() --> 697 update = self.process_update( 698 path, 699 process, 700 state, 701 process_interval 702 ) 704 # update front, to be applied at its projected time 705 self.front[path]['time'] = future File ~/code/process-bigraph/process_bigraph/composite.py:636, in Composite.process_update(self, path, process, states, interval, ports_key) 614 """Start generating a process's update. 615 616 This function is similar to :py:meth:`_invoke_process` except in (...) 631 ``store``. 632 """ 634 states = strip_schema_keys(states) --> 636 update = process['instance'].invoke(states, interval) 638 def defer_project(update, args): 639 schema, state, path = args File ~/code/process-bigraph/process_bigraph/composite.py:303, in Process.invoke(self, state, interval) 302 def invoke(self, state, interval): --> 303 update = self.update(state, interval) 304 sync = SyncUpdate(update) 305 return sync File ~/code/spatio-flux/spatio_flux/processes/particles.py:197, in Particles.update(self, state, interval) 194 x, y = get_bin_position(new_position, self.config['n_bins'], self.env_size) 196 # Update local environment values for each particle --> 197 updated_particle['local'] = get_local_field_values(fields, column=x, row=y) 199 # Apply exchanges and reset 200 exchange = particle['exchange'] File ~/code/spatio-flux/spatio_flux/processes/particles.py:42, in get_local_field_values(fields, column, row) 40 local_values = {} 41 for mol_id, field in fields.items(): ---> 42 local_values[mol_id] = field[column, row] 44 return local_values KeyError: (1, 30)
In [ ]:
v6x.plot_particles_snapshots(skip_frames=3)
Particle-COMETS¶
In [46]:
dfba_config = {
"model_file": "textbook",
"kinetic_params": {
"glucose": (0.5, 1),
"acetate": (0.5, 2)},
"substrate_update_reactions": {
"glucose": "EX_glc__D_e",
"acetate": "EX_ac_e"},
"biomass_identifier": "biomass",
"bounds": {
"EX_o2_e": {"lower": -2, "upper": None},
"ATPM": {"lower": 1, "upper": 1}
}}
In [47]:
bounds = (10.0, 20.0) # Bounds of the environment
n_bins = (4, 8) # Number of bins in the x and y directions
v7 = SpatioFluxVivarium()
v7.add_process(
name='particle_movement',
process_id='Particles',
config={
'n_bins': n_bins,
'bounds': bounds,
'diffusion_rate': 0.1,
'advection_rate': (0, -0.1),
'add_probability': 0.1,
'boundary_to_add': ['top']})
v7.connect_process(
name='particle_movement',
inputs={'fields': ['fields'],
'particles': ['particles']},
outputs={'fields': ['fields'],
'particles': ['particles']})
# add fields and diffusion process
mol_ids = ['glucose', 'acetate', 'biomass']
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {'biomass': (0, -0.1)}
# initialize the molecular fields
max_glc = 10
glc_field = np.random.rand(n_bins[0], n_bins[1]) * max_glc
acetate_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field = np.zeros((n_bins[0], n_bins[1]))
biomass_field[0:int(1*n_bins[0]/5), int(2*n_bins[1]/5):int(3*n_bins[1]/5)] = 0.1 # place some biomass
v7.add_object(name='glucose', path=['fields'], value=glc_field)
v7.add_object(name='biomass', path=['fields'], value=biomass_field)
v7.add_object(name='acetate', path=['fields'], value=acetate_field)
# add a spatial diffusion advection process
v7.add_process(
name='diffusion_advection',
process_id='DiffusionAdvection',
config={
'n_bins': n_bins,
'bounds': bounds,
'default_diffusion_rate': diffusion_rate,
'default_diffusion_dt': diffusion_dt,
'advection_coeffs': advection_coeffs},
inputs={'fields': ['fields']},
outputs={'fields': ['fields']})
# add a dynamic FBA process at every location
for i in range(n_bins[0]):
for j in range(n_bins[1]):
dfba_name = f"dFBA[{i},{j}]"
v7.add_process(
name=dfba_name,
process_id="DynamicFBA",
config=dfba_config)
v7.connect_process(
name=dfba_name,
inputs={
"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids}},
outputs={
"substrates": {
mol_id: ['fields', mol_id, i, j]
for mol_id in mol_ids}})
v7.initialize_process(
path='particle_movement',
config={'n_particles': 2})
v7.add_emitter()
v7.diagram(dpi='70',
remove_nodes=[
f"/dFBA[{i},{j}]"
for i in range(n_bins[0]-1)
for j in range(n_bins[1])]
)
Out[47]:
In [48]:
v7.run(60)
In [49]:
v7.plot_particles_snapshots(skip_frames=3)
Saving GIF to species_distribution_with_particles.gif
